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A Note on the Numerical Integration of Differential 

Equations 1 

By W.E.Milne 2 

An integration method for ordinary differential equations is developed, in which the 
approximation formulae contain derivatives of higher order than those contained in the 
differential equation itself. The method is particularly useful for linear differential equa- 
tions. Numerical examples are given for Bessel's* differential equation. 



I. Introduction 

The object of this note is to present a method for 
the numerical integration of ordinary differential 
equations that appears to possess rather outstand- 
ing advantages when applied to certain types of 
equations. The equations to which the method 
most readily applies are those for which it is pos- 
sible to obtain, in comparatively simple form, 
expressions for two additional derivatives. That 
is, for an equation of n-th. order 

we obtain by differentiation expressions for y {n+1) 
and 2/ (n+2) . If these expressions are not so 
involved as to make the labor of substitution 
prohibitive, then the method here proposed is 
applicable. 

The advantages claimed for the process are: 

(1) The start of the integration is accomplished 
by the same formulas that are used in the regular 
routine of the process, so that no special formulas 
or procedures are required in order to get the 
computation underway. 

(2) Each step of the integration makes use of 
only two lines of the computation, whereas a 
method employing differences and having a com- 
parable degree of accuracy would require five lines 
in the computation. 

(3) A change in the length of step for step-by- 
step integration is often necessary as the integra- 

1 The preparation of this paper was sponsored (in part) by the Office of 
Naval Research. 

2 Mathematics Department, Oregon State College. 



tion proceeds. Such a change can be more readily 
made in this process than where five-line formulas 
of integration are employed . 

(4) The coefficients occurring in the formulas are 
simpler than those in comparable five-line quad- 
rature formulas, so that the machine calculation is 
not at all complicated. 

The most obvious disadvantage of the process 
is that it requires the calculation of two additional 
derivatives at each step, and the labor of substi- 
tution in certain instances may be excessive. In 
such cases this method is not recommended. On 
the other hand, for equations of simple analytical 
form, and particularly for linear differential equa- 
tions, it should prove valuable. 

II. Derivation of Formulas 

Let x and x u where Xi — x =h, be two values 
of the independent variable x, and let y 0} Vo, etc.. yiy[ 
etc., be the corresponding values of y, y' , . . . As- 
suming that y has a continuous derivative of order 
7 we may express y, y f , y" , and y nr by Taylor's 
series with remainder term, as follows: 



hhi hh/ i3) h 4 v i4) // 5 7/ (5) /> 6 ?/ (6) 

n, — lt i hn f I y ° I- y ° 1 y ° -r- y ° -4- y ° 
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/,5, y (5) J, 6^.(6) 

y2=yo+2%o+-o,--h^T g -+ 7i 



2! 
32A 8 yff» , 64A 6 yg 



5! 



6! 
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4! 



(4) 



(5) 



From the four equations (1, 2, 3, and 4), the three 
quantities y ( 4) , y ( 5) , and t/o 6) may be eliminated, and 
the resulting equation can be rearranged so as to 
give 

h h 2 h 3 

yi-yo=2 (2/1+2/0) — jq (y"i-yl) + ^ toA' +&') + 



Rfi 



(6) 



It may be shown by a separate investigation 3 
that 

g -A y y (7) (s ) 
6 100,800 ' 

in which x <s<^i. 

In a similar manner from eq 1 to 5 we may de- 
rive 

y 2 -2y l +y =7h(y[-y' )-Zh 2 (y:+y'Z) + 



h3 m '" k »'\ ■ 210ft 7 y (7) W 
12 < U * ~ 5 ^ > + -T6p00-' 



(7) 



with 8 in the interval x <s<x 2 . These are the 
required formulas. 

III. Application to Equations of First Order 

Let the given differential equation be 

v'=Mv)- (8) 

Differentiation gives 

V"=Mx,y)+M*,V)v', 0) 

and 

y'"=.fxx(x,y)+2f xy (x,y)y'+f„(x,y)y ,2 +j u (x,y)y". 

(10) 

Let the given initial values be x , y . Then 
from eq 8, 9, and 10 in succession we obtain 
Vo, y'o, y'o" , giving the first line of the computation 



x T/o 



Vo 



Vo 



in 

y* . 



3 W. E. Milne, The remainder in linear methods of approximation, J. 
Research NBS 43, 501 (1949) RP2042. 



To proceed we assume a trial value for y x . A 
fairly good trial value is provided by 

Next with x h and the trial value y u we obtain 
trial values for y[, y" , y" from eq 8, 9, and 10 
and have the trial line for X\ : 



Xl 



Vi 



Vi 



Vi 



Vi 



Now using eq 6 we secure an improved value of 
y 1} compute y[, y[\ y[" , from eq 8, 9, and 10, 
recalculate y 1 by 2, (eq 6), and repeat this sequence 
of steps until no change occurs in the value ob- 
tained for y x . This is taken as correct, and we 
have two lines of the computation: 



2*0 


yo 


yo 


yo 


yo 


Xi 


yi 


y'i 


v" 


y'i 



We are now in a position to use formula 7 in 
order to calculate a trial value for y 2 . Then trial 
values of y' 2) y^ ', y*' are obtained from eq 8, 9, 
and 10, and we leave the trial line: 



x 2 



2/2 



2/2 



2/2 



fit 

y 2 . 



Formula 6 (with subscripts advanced by 1) gives 
an improved value for y 2 . If the " improved" 
value of y 2 is different from the "trial" value of 
y 2 , it will be necessary to recompute y f 2 , y 2 , y 2 ' 
and apply eq 6 again, repeating these steps until 
no further change occurs. At this stage we have 
completed three lines of the computation: 



Xq 


yo 


yo 


yo 


y , 


Xi 


yi 


y'i 


Vi 


ft r 
Vl , 


x 2 


2/2 


2/2 


v% 


2/2". 



Subsequent lines of the calculation are obtained 
by exactly the same steps as were used to get y 2 . 

IV. Discussion of the Process 

This completes the description of the process 
for the case of equations of the first order. Some 
comments of a practical nature are, however, 
pertinent. 

(1) Obviously, the error in y at any step due to 
the use of the approximate formula 6 is bounded 
by the quantity 

K ~mm' (11) 
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in which .M=max|?/ (7) | in the interval covered by 
the step. Supposing for the moment that M is 
a known constant, we see that the above equation 
connects the magnitude of the error E with the 
length of the step-interval h. 

If, for example, the problem in question requires 
that y be obtained accurately to a specified num- 
ber of decimal places, eq 11 enables us to select an 
appropriate value of h that will secure this ac- 
curacy. On the other hand, if we have already 
decided on the value of h, eq 11 will tell us how 
many decimals in the result may be regarded as 
correct. 

(2) Actually, M is rarely constant from step 
to step, and moreover the value of M is unknown 
since it is ordinarily utterly impractical to calcu- 
late the value of the seventh derivative of y at 
each step. However, a crude estimate (which 
actually proves to be sufficiently satisfactory) 
may be made as follows: Assume that the calcula- 
tion has been performed correctly to the nth step 
so that we have the correct values of the line 



Vn 



Vn 



Vn 



Vn 



a trial value y n +\ is obtained by eq 7. A final 
value of y n +i is obtained by eq 6, repeated if 
necessary. Now the error of y n +\ is 

2iohy 7) (s) 

100,800 ' 
whereas that of y n +i is 

fe¥ 7) (Q , „ 



100,800 



e being the error produced in y n +\ by the errors in 
y'n+i, y'n\-u y'n'-lh these latter being due to the 
error of y n +i. In actual practice, if the process is 
rapidly enough convergent for practical use, the 
error e must be much smaller than the error in 
2/n+i. Hence, we may neglect e. Then 



VnH Vn+l z 



2iohy 7) (s) . >y 7 >(o 

100,800 *" 100,800* 



If we ignore the fact that y i7) (s) and y i7) (s') are 
not exactly the same, we may add the terms on 
the right, divide by 211, and obtain 



Error of y n 



h 7 y (7) (s)_y n+l —y n+1 



Although the foregoing reasoning seems very 
crude, the final formula 

Error of y„ +1 = ym+1 ~f" +1 , 

proves to be not only simple in application but 
actually quite reliable in determining how many 
significant figures of the result can be trusted. 

(3). At each step of the computation the quan- 
tity c n =y n —y n should be recorded in a separate 
column. This column of c's is used for several 
purposes. 

(a) As long as the c's vary regularly and have 
significant figures only in the last two places 
retained, we proceed with the computation in 
reasonable confidence that all is well. 

(1)) A sudden fluctuation in the c's suggests 
that a computational mistake has been made, and 
the lines involved should be rechecked. 

(c) If the c's increase to the point of affecting 
the last three places retained, then either the 
interval, h, should be shortened, or one less place 
should be retained. 

(d) The necessity for recomputing y' n +u Vn+iy 
y'nli can frequently be obviated by estimating 
c n+l from the known values c n _ 2 , c n _i, c n and add- 
ing it to the trial value y n+1 before computing 
y n +i, pn+i, Vn+i- If the estimate is sufficiently 
accurate, no recomputation is required. 

(4). The foregoing discussion applies only after 
the computation is under way and does not give 
any clue to the accuracy of y { . We would, of 
course, like to decide on the value of h and 
on the number of decimal places before starting 
the computation. This requires that we calculate 
y ( 7) from the differential equation and the initial 
values. 

V. Equations of Higher Order 

The modifications required to apply the process 
to equations of second or higher order are slight. 
In the case of an equation of second order, for 
example, the routine (after the start has been 
made) is as follows: 

Predict y' 2 by eq 7 modified as follows: 



100,800 



211 



y2-2y[ + y f =7h(y f 1 f -y , f )- 
§ (ll<-52/n. 



-3W+2/D + 
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Predict y 2 by eq 7. Calculate y 2 f , y' 2 " , y^ from 
the differential equations and the equations 
obtained by differentiation. 

Correct y' 2 by 

h h 2 h 3 

Correct y 2 by 



VI. Illustrative Examples 

The foregoing method is applied to the second 
order differential equation 

xy"+y'+xy=0, 

with conditions y=l, y /=z at x=0. 

Example 1, below, gives the solution by another 
method using h=0.1. Example 2 uses the present 



method with h =0.5. A comparison of the com- 
putation time is given. Example 3 uses the pre- 
sent method with h =0.1. 

It appears that in example 2 the error of y is 
about 2 in sixth place, whereas in example 3 the 
error is occasionally 1 in the tenth place. 

A comparison of example 1 and example 2 
shows that the new method obtained the value of 
Jo (3) in only six steps (and more accurately) 
than the simple method based on Simpson's 
Rule could secure in 30 steps. Although the 
labor of substitution per step is much greater for 
the new method, the reduction in number of re- 
quired steps more than offsets this extra work, as 
is shown by a comparison of required times. 

Example 3 provides further evidence of the 
power of the new method. Anyone with experi- 
ence in numerical solution of differential equations 
will recognize that to solve Bessel's equation to 
10 decimal places in the neighborhood of the orign 
with step-intervals of length 0.1 requires a pretty 
accurate method of solution. 



540 



Journal of Research 



Example 1. 






dty <hi 
Differential equation: x -j-^-\-^--{-xy=0 
dx 2 dx 

Computation formulas 

Predictor : y' n+l = 2/^_ 3 +/i[4y^_ 1 +8/3 «2j£_ x ] . 

Corrector: yi +J =jC-i+A[2l£+l/3 &y' n ], 

1 
Derivatives: y" = y' — y. 

Time: 2H hr; A=0.1 











l A W 


n w 






True values 




y 


V 


y 






y 


v' 


0.2 
.3 
.4 
.5 
.6 


»0. 990025 
a. 977626 
ft. 960398 
a. 938470 
. 912005 


a -0.099501 
»-. 148319 
»-. 196027 
a-. 242268 
-.286702 (0) 


-0. 492520 
-.483229 
-.470330 
-.453934 
-.434168 


0. 000370 
. 000489 
. 000603 
. 000714 







2 


0. 912005 




0.001203 
.001166 
.001122 
. 001065 








-0.286701 


.7 


.881201 


-.328995 


-.411208 


.000815 


.001004 








.881201 


-.328996 


.8 


. 846288 


-.368843 (2) 


-.385234 


.000914 


. 000930 





1 


. 846287 


-.368842 


.9 


. 807524 


-.405949 (8) 


-.356470 


. 001001 


. 000853 





1 


. 807524 


-.405950 


1.0 


. 765198 


-.440052 


-.325146 


. 001085 


. 000764 








. 765198 


-.440051 


1.1 


. 719622 


-.470902 (1) 


-.291529 


.001154 


. 000674 





1 


. 719622 


-.470002 


1.2 


b. 671133 (5) 


-.498290 


-.255891 


.001219 


. 000574 


2 





.671133 


-.498289 


1.3 


. 620086 


-.522023 (1) 


-. 218530 


.001269 


. 000473 





2 


. 620086 


-.522023 


1.4 


.566855 (7) 


-. 541949 


-.179749 


. 001313 


. 000366 


2 





. 566855 


-.541948 


1.5 


. 511828 


-. 557936 


-. 139871 


. 001342 


. 000259 








.511828 


-.557937 


1.6 


.455402 (4) 


-. 569897 


-.099216 


. 001364 


. 000146 


2 





. 455402 


-.569896 


1.7 


. 397985 


-.577765 (4) 


-.058123 


. 001372 


+.0000:? 7 





1 


. 397985 


-.577765 


1.8 


.339986 (7) 


-.581518 (9) 


-.016920 


. 001371 


-.000077 


1 


-1 


. 339986 


-.581517 


1.9 


.281819 (20) 


-.581157 (6) 


+. 024053 


. 001357 


-.000184 


1 


1 


. 281819 


-.581157 


2.0 


.223890 (2) 


-.576726 (7) 


. 064473 


. 001334 


-.000295 


2 


-1 


. 223891 


-.576725 


2.1 


.166607 (8) 


-.568292 (1) 


. 104008 


.001298 


-.000398 


1 


+1 


. 166607 


-.568292 


2.2 


.110361 (3) 


-.555964 (5) 


. 142350 


. 001255 


-.000502 


2 


_! 


. 110362 


-.555963 


2.3 


.055540 (1) 


-.539872 


. 179187 


.001198 


-.000595 


1 





. 055540 


-.539873 


2.4 


.002506 (8) 


-.520186 (8) 


. 214238 


.001136 


-.000689 


2 


-2 


+. 002508 


-.520185 


2.5 


-.048384 (3) 


-.497093 


. 247221 


. 001060 


-.000771 


+1 





-.048384 


-.497094 


2.G 


-.096807 (5) 


-.470819 (20) 


. 277891 


. 000982 


-.000852 


+2 


-1 


-.096805 


-.470818 


2.7 


-.142450 (49) 


-.441600 


. 306006 


. 000890 


-.000919 


+1 


0. 


-. 142449 


-.441601 


2.8 


-.185038 (6) 


-.409710 (1) 


. 331363 


. 000798 


-.000984 


+2 


-1 


-. 185036 


-.409709 


2.9 


-. 224312 


-.375426 


. 353769 


. 000694 


-.001034 








-.224312 


-.375427 


3.0 


-.260054 (3) 


-.339060 (1) 


. 373074 






+1 


-1 


-.260052 


-.339059 





» Starting values given. 

b In the column for y and y' appear the corrected values. 



The digits of the predicted values when different from the corrected are shown in parentheses. 
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Example 2. 



d%y dy 
Differential equation: x — ^-|- — -{-:ry=0 

Computation formulas 



ft 3 



Predictor: fi+i-fci-iri-i+Wfri'-fiLi) ~3^ 2 ^n '+^i) +^ {pff~ti5-i> 



y n +i = ' 2 yn-yn-i^h{y f rl -y f n _ l )-ZhHy' y l J ty , n-i) + — (ll^'-Sy^O 

/i h 2 /i 3 

Corrector: K l+ i=^ + v W+i+C)-^ (vi+l-»n')+^ M+i+?» v ) 






^ ft 2 ft 3 



1 , 

Derivatives: 2/ =-- y 



2 „ , 1 



120 



2/ IV = y'"-y' 

x 

Time: 1H hr; /i=0.5 



a; 


2/ 


y' 


y" 


y'" 


ylV 


c 


c' 


True values 


y 


y' 




0.5 

1.0 

1.5 

2.0 

2.5 

3.0 


•1. 

0. 938470 
b . 765195 (9) 
. 511826 (31) 
. 223889 (2) 
-.048382 (6) 
-.260053 (47) 


a 

a— 0. 242268 
-.440047 (63) 
-. 557934 (5) 
-.576721 (18) 
-.497090 (1) 
-.339057(48) 


-0 5 
-.453934 
-.325148 
-. 139870 
-h 064472 
. 247218 
.373066 




0. 181064 
. 325148 
. 403210 
. 400304 
. 318668 
.177021 


0.375 

. 336622 

. 229798 

. 077382 

-.088207 

-.231948 

-.324043 


4 

5 

-7 

-4 

6 


-16 
-1 

+3 
-6 
+3 


0. 765198 

.511828 

. 223891 

-. 048384 

-.260052 






-0. 440051 
-. 557937 
-.576725 
-.497094 
-.339059 



a Starting values given. 

b In the column for y and y' appear the corrected values. The digits of the predicted values when different from the corrected are shown in parentheses. 



Example 3. 



d2 y . d y r 

Differential equation: x ~r^~r^r ~r xy=0 m 

Computation formulas: (Same as Example 2.). 
Computation, ft=0.1. 



















True values 




y 


y' 


y" 


y'" 


ylV 


c 


c' 


y 


y' 


0.5 


H 


»0 


-0.5 





0.375 












1 


a 99750 15621 


a-0. 04993 75260 


—.49812 63017 


0. 03744 79394 


.37343 86390 












.2 


. 99002 49723 


-.09950 08326(7) 


-. 49252 08093 


. 07458 40641 


. 36876 81738 


-0 


-1 


0.2 


0. 99002 49722 


-0. 09950 08326 


.3 


. 97762 62466 


-. 14831 88162 (6) 


-.48323 01926 


.11109 92782 


. 36102 95182 





-4 


.3 


. 97762 62465 


-. 14S31 88163 


.4 


. 96039 82267 


-.19602 65779 (81) 


-.47033 17820 


. 14668 99212 


. 35029 02625 





-2 


.4 


. 96039 82267 


-.19602 65780 


.5 


b. 93846 98073 (2) 


-.24226 84576 (9) 


-. 45393 28921 


. 18106 04114 


. 33664 42541 


-1 


-3 


.5 


. 93846 98072 


-.24226 84577 


.6 


.91200 48636(5) 


-.28670 09880 (1) 


-.43416 98836 


. 21392 58273 


. 32021 07071 


-1 


-1 


.6 


. 91200 48635 


-.28670 09881 


.7 


. 88120 08887 


-.32899 57415 (6) 


-.41120 69723 


.24501 43928 


. 30113 31217 





-1 


.7 


. 88120 08886 


-.32899 57415 


.8 


. 84628 73528 (7) 


-.36884 20461 (2) 


-.38523 47952 


. 27406 98431 


. 27957 79988 


-1 


-1 


.8 


. 84628 73528 


-.36884 20461 


.9 


. 80752 37982 (0) 


-.40594 95461 (2) 


-.35646 87470 


. 30085 36525 


. 25573 33411 


-2 


-1 


.9 


. 80752 37981 


-.40594 95461 


1.0 


. 76519 76866 (0) 


-.44005 05858 (9) 


-.32514 71008 


. 32514 71008 


. 22980 69700 


-6 


~ '1 


1.0 


. 76519 76866 


-. 44005 05857 



» Starting values given. 

b In the column for y and y' appear the corrected values. The digits of the predicted values when different from the corrected are shown in parentheses. 
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